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ABSTRACT 

The discovery of multi-planet extrasolar systems has kindled interest in using their orbital evolution 
as a probe of planet formation. Accurate descriptions of planetary orbits identify systems which could 
hide additional planets or be in a special dynamical state, and inform targeted follow-up observations. 
We combine published radial velocity data with Markov Chain Monte Carlo analyses in order to obtain 
an ensemble of masses, semimajor axes, eccentricities and orbital angles for each of 5 dynamically active 
multi-planet systems: HD 11964, HD 38529, HD 108874, HD 168443, and HD 190360. We dynamically 
evolve these systems using 52,000 long-term N-body integrations that sample the full range of possible 
line-of-sight and relative inclinations, and we report on the system stability, secular evolution and the 
extent of the resonant interactions. We find that planetary orbits in hierarchical systems exhibit 
complex dynamics and can become highly eccentric and maybe significantly inclined. Additionally 
we incorporate the effects of general relativity in the long-term simulations and demonstrate that can 
qualitatively affect the dynamics of some systems with high relative inclinations. The simulations 
quantify the likelihood of different dynamical regimes for each system and highlight the dangers of 
restricting simulation phase space to a single set of initial conditions or coplanar orbits. 
Subject headings: celestial mechanics — planetary systems: formation — methods: N-body simula- 
tions, statistical 



1. INTRODUCTION 

Currently, over 30 multi-planet exosystems are each 
known to include 2-5 known planets. The orbital 
architectures and formation scenarios for these sys- 
tems have become subjects of numerous investiga- 
tions. T he recent discovery of a fifth planet around 
55 Cnc (|Fischer et all [20081) has prompted a flurry 
of follow-up studies dGavon et al.1 [20081 IRavmond et all 
I2008at Ui et al.1 |2QQ9j that aim to better describe 
the dynamics and evolution of that system. De- 
spite sparse data, the dir ectly-imaged triple system HR 
8799 dMarois et al.l I2008D has bee n sub ject to intense 
scrut i ny (iFabrvckv fe Murrav-Clavll2008t Iciose fc Males] 
; iFukagawa et al.ll2009l; IGozdziewski fe Migaszew ski 



lLafreniere et all 120091: iReidemeister et all "12009; 



Sudol et al.ll2009f) . These explorations demonstrate that 



major open questions regarding multi-planet systems re- 
main, including: What is the timescale for instability in 
these systems? Could they admit additional, currently 
undetectable planets? How tightl y "packed" are such 
systems (e.g. IRavmond et al.ll2009D ? 

These questions can be addressed by investigating the 
orbital evolution of individual planetary systems. How 
do the orbital eccentricities and semimajor axes of plan- 
ets change with time? How does additional physics (e.g. 
tidal forces and general relativity, henceforth referred to 
as GR) affect the orbits of short-period planets and in- 
directly other planets in the system? Theoretical inves- 
tigations addressing these questions should account for 
the uncertainties in orbital parameters obtained from ob- 
servations. In order to study dynamics, we need to es- 
tablish initial conditions. Unfortunately, the accuracy 
of such initial conditions are limited due to 
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ment uncertainties and degeneracies inherent to the ra- 
dial velocity (RV) exoplanet discovery technique. Oc- 
casionally, the best-fit RV data yields parameters which 
indicate that the timescale for such a system to undergo 
instability is much less than the age of the system (e.g. , 
in less than 10 5 yr, for HP 8 2943; iMavor et all 120041 : 
IFabrvckv fc Murr av-Clavl 120081 ) . These short timescales 
(relative to system lifetimes) suggest that the best-fit or- 
bital model is unlikely and motivate investigations to find 
the plausible and stable solutions. Further, both i re i and 
ilos for planets in most multi-planet systems have not 
yet been measured. Noteworthy exceptions include a re- 
cent esti mate of the relative in clination (i re i) of GJ 876 
planets (jBean fc Scifahrt 200^) and recent estimates for 
the li ne-of-sight inclinations (ilos) of the planets in HR 
8799 (jLafreniere et al.ll2009D . 

In an effort to remedy some of these issues and bet- 
ter describe system properties, investigators have been 
developing techniques to mo del RV data and to de- 
scri be their orbital evolution . IGozdziew ski et all ()2005f ) 
and lGozdziewski"et al. (2008) utilize stability constraints 
and optimization techniques partly derived from the 
MEGNO (Mean Exponential Growth Fac tor for Nearby 
Orbits) method (|Cincotta fc Simol 120001 ) . The advan- 
tages to MEGNO include simultaneous multi-planet fit- 
ting and efficient identification of quasi-periodic or ir- 
regular (chaotic) motion. A disadvantage is that the 
arbitrary choice of the "penalty" parameter determin- 
ing this timescale prevents a rigorous i n terpretation of 
the re sults. Some authors (|Fordl 120051 120061: iGregorvl 
I2007al|bl ) use Markov Chain Monte Carlo (MCMC) tech- 
niques to generate a sample of initial conditions from the 
posterior probability distribution. The strength of these 
techniques relies on rigorous Bayesian calculation and in- 
terpretations, and a weakness is that stability must be 
individually tested in each of a large number of mod- 
els. For datasets that result in many unstable models, 
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this procedure may result in a time-consuming and per- 
haps inefficient co mputational effort. Here, we follow 
the methodology of lVeras fe Fordl (|2009t ). which relies on 
MCMC-based simulations to derive ensembles of initial 
conditions. 

We assume that the motion is described by a sum of 
Keplerian (i.e. non-interacting) orbits. Each ensemble of 
initial conditions consists of a list of masses, semimajor 
axes, eccentricities, arguments of pericenter, inclinations 
and nodes, assuming an initial mean anomaly for each 
planet. We then assign line-of-sight and relative inclina- 
tions to planets with these ensembles of orbital elements 
in order to sample from the entire phase space of possible 
initial conditions. 

When planetary orbits are integrated forward in time, 
the system may exhibit a variety of evolutionary paths. 
Unrestricted planetary inclinations admit a wide vari- 
ety of phase space regimes, often including those where 
apsidal and resonant angles circulate, librate and ex- 
hibit chaotic motions (fMichtchcn ko et aDl2006bh . In this 
work, because we consider systems containing two mas- 
sive exoplanets with a wide range of eccentricities and 
inclinations, we must appeal to numerical integrations. 

Here, we investigate the dynamical evolution of five 
two-planet systems which contain sufficiently numerous 
and accurate radial velocity orbital data suitable for a 
wide-ranging study. In $31 we show that the stability 
and dynamical properties of these systems can be signif- 
icantly influenced by the initial values of i re i and ilos- 
Four of these systems are "hierarchical", meaning they 
contain large ratios of orbital distances, which often in- 
clude a close-in planet whose evolution could be affected 
by the general relativistic precession of the pericenter. 
We aim to characterize the eccentricity variation, stabil- 
ity, secular effects and resonant signatures of planets in 
these systems as a function of i re i and ilos, while tak- 
ing into account the uncertainties of the measured orbital 
parameters for each planet. We combine investigations 
of the hierarchical HD 12661 with the system studied 
here to draw inferences for effective methods of future 
dynamical investigations of hierarchical systems, and to 
determine their general evolutionary trends (such as the 
propensity, or lack t hereof, for planets to period ically at- 
tain circular orbits: [Barnes fe Greenberg1l2006al lb1). 

We describe our methodology in and present the 
results of the N-body simulations in |j3j Tables [2"][T4l 
present summary statistics, two or three tables per sys- 
tem, arranged by row according to the binned relative 
inclination between both planets; the accompanying fig- 
ures help the reader visualize the data. We discuss the 
results in g]and conclude in $5j 

2. METHODS 

The setup o f our simulations closely follows that of 
iVeras fe Fordl (|2009l ). First, we generate a sample of ini- 
tial conditions that will serve as the basis for numerical 
integrations. 

Using RV data from lWright et all (|2008f ). we utilize the 
MCMC an alysis methods and prior distributions from 
iFordl (|2006|) . In particular, we calculate 5 Markov chains 
each containing 10 5 or 10 6 states per system. Each state 
includes the orbital period (P), velocity amplitude (K), 
eccentricity (e), argument of pericenter measured from 
the plane of the sky (w), and mean anomaly at a given 



epoch (u) for each planet specified by subscripts. We 
randomly select 13 sets of 500 states, where each set is re- 
stricted to a particular range of i re i. We sampled from a 
uniform distribution for the longitude of ascending nodes, 
Q, and from an isotropic distribution for the ilos of each 
planet. Then we used rejection sampling to choose 500 
isotropically distributed pairs of ilos values in each of 
13 i re i bins. The planet masses, m, and semimajor axes, 
a, are obtained from each set of (P, K 1 e, u), i, Q, u) values 
from relations deri ved with a Jacobi coordinate system 
(|Lee fc Peald 12003). The 13 i re i bins are split into one 
coplanar bin and 12 bins which contain systems with rel- 
ative inclinations in 15° intervals from 0° to 180°. We 
used stratified sampling to obtain sufficiently large sam- 
ples at both low and high inclinations. 

Second, we perform integrations to test each set of 
initial conditions for long-term stability. We integrated 
each set of initial conditio ns with the hyb rid symplec- 
tic integrator of Mercury (|Chambersl Il999f ) for 1 Myr, 
and obtained snapshots of the system states every 10 3 
yr. Inspection of a representative sample of stable con- 
figurations reveals that 1 Myr typically exceeds secular 
timescales. 

We assume that the actual systems will not undergo in- 
stability immediately (< 10 6 yr), and generate a subsam- 
ple of initial conditions that do not show indications of in- 
stability during our simulations. We classified systems as 
"unstable" if, for either planet, (a max — £imin)/»o > to-o, 
where a max , a m j n and ao represent the maximum, mini- 
mum and initial values of the semimajor axis, and r = 
0.3. Additionally we consider systems which have under- 
gone a close encounter (where two bodies come within a 
Hill radii of one another) or collision to be unstable. Our 
criterion may not identify a small fraction of some sys- 
tems that will manifest instability on longer timescales, 
but importantly it also avoids miscategoriz ing a system 
exhib iting bounded chaos as unstable (e.g., Gozdziewski 
l200l . 

In addition to testing for stability, we compute the ex- 
tent to which planets' eccentricities and mutual relative 
inclination vary over the course of the simulations and 
for what fraction of the planets' eccentricity and incli- 
nation exceed their initial values. We also determine 
how closely the planetary orbits approach circularity by 
computing K b = e^ min /e b ,max, where e b , min and e b , m ax 
represent the minimum and maximum eccentricities at- 
tained throughout a simulation for planet b. A similar 
definition holds for planet c. If k ~ 0, a planet's eccen- 
tricity will periodically approach zero. 

For two planets on nearly planar orbits, the differ- 
ence of the longi tude of periapses, Am, has dynam- 



ical si gnificance (IChiang et al.l 120011: IChiang fc Murray 
20021 : iBeauge et all 120031: Ui et all 120031: IZhou fc Sun 



27)031 : In amou nil l2005l : iMigaszewski fc Go 

!dziewskl2009) 
In the planar case, when this difference is 0° 
(180°), the planetary orbits are said to be "aligned" 
( "anti- aligned"). An "Apsidal Cor o tation Resonance" 
(ACR : IBeauge fc Michtchenkol 120031 iFerraz-Mello et al l 



2001 iLed 12004: iKlev et all 120051: Beauge et al.l 120061: 



Mich tchenko et al. |2006eJ: Voyatzis & Hadiidcmetriou 
20061 : iSandor et al.l l2007t iMichtchenko et al.1 l2008al lbl) 
refers to system evolutions that result in Atu librat- 
ing about 0° or 180° and can result from either purely 
secular evolution or secular plus resonant interactions. 
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However, for nonzero inclinations, the angle between the 
pericenter directions deviates from Aw in the planar 
case. The most natural plane on which to measure three- 
dimensional dynamics is the invariable plane. Therefore, 
for each snapshot in time and for each N-body integra- 
tion, we measure Aw', the difference of the longitude of 
pericenters after projecting onto the invariable plane. A 
time series of these measurements of Aw' over a period 
of time, ioj sometimes reveals either "circulation," or "li- 
bration" about a particular value, the "libration center," 
^o- Strictly, an angle librates if there is a range of val- 
ues that the angle avoids; Iq is then the middle value 
of the range not avoided by the angle. Unless the an- 
gle is continuously sampled, there will always be a finite 
range avoided by the angle. Therefore, because of finite 
sampling from integration output, we consider Aw' to li- 
brate if every value of Aw' over a time to avoids a range 
larger than 90°. The location of this range depends on 
the value of lo considered. 

When testing for libration, we consider a putative value 
of Iq and determine if every value of Aw' is contained in 
the range from lo — 135° to Iq + 135°. If an angle librates, 
inv estigators often rep ort a libration amplitude about 
lo- IVeras fc For d (2009) utilize two alternative variation 
measures, the root mean square deviation and mean ab- 
solute deviation, which are more robust than the am- 
plitude given finite sampling and additional short-term 
perturbations. We calculate the RMS deviation of Aw' 
about each of l = (0°, 90°, 180°, 270°). Our choice of l 
includes each of th e common libration cente r s. Analyt- 
ical c onsiderations ([Murray fe Dermottlll999l; iMorbidellil 
I2002D demonstrate that l = 0° and l = 180° ("symmet- 
ric solutions" ) are the only possiblities for purely secular 
evolution, and numerous numerical studies of aligned and 
anti-aligned configurations bear this trend out. However, 
lo may take on other values ("a s ymme t ric solutions" 



to may take on otner values ( asymmetric solutions ; 
Beaugel 119941: IWinter fc Murray! 1l997t Pancart et al.1 



2002HBeauge et all 120031: iMurrav-Clav fc Chiang! 12005^ 
if a system is near a mean motion resonance. Of the sys- 
tems considered here, only HD 108874 is near a MMR. 
Assuming the system is not undergoing asymmetric libra- 
tion and considering a few alternative values of lo (e.g., 
90° and 270°) allows us to estimate the frequency of sys- 
tems being misclassified. 

3. RESULTS 

3.1. Common Attributes 

Wri ght et al.l (|2008f ) provide best-fit orbital solutions 
to radial velocity observations of multiplanet systems. 
Our calculations are based on the posterior probabil- 
ity distribution for orbital parameters described in Sec- 
tion [21 Here, we give an approximate summary of 
system properties. In all systems except HD 190360, 
planet b is the inner planet. The ratio of semima- 
jor axes nearly exceeds one order of magnitude in HD 
11964 (a b ,a c « 0.2,3.2 AU), HD 38529 (a b ,a c 
0.1,3.7 AU), HD 168443 {a b ,a c « 0.3,2.9 AU), and HD 
190360 (a Cl a b ~ 0.1,3.9 AU), and approximately cor- 
responds to the 4:1 Mean Motion Resonance (MMR) 
in HD 108874 (a b ,a c « 1.1,2.7 AU). Assuming copla- 
narity, The outer planet is more massive in HD 11964 
(M b sin i hLO s, M c sin i c , LO s « 0.1, 0.7M Jup ), HD 38529 
(M b sin i b , LO s,M c sin i c . LO s « 0.8,13M Jup ), HD 168443 



(M b smi b , LO s, M c sin i c ,Los ~ 8, 18Af Jup ), and HD 
190360 (M c smi CtLO s,M b 8mi b>LO s « 0.06, 1.5M Jup ), 
whereas the planet masses are comparable in HD 108874 
(M b sin i b ,LOS, M c sin i c ,LOS ~ 1.4, 1.0M Jup ). In the 
limit that a planet's orbit is seen face-on, the planet's 
mass diverges and the corresponding system becomes un- 
stable. However, most orbital orientations with respect 
to the Earth yield actual planet masses less than twice 
the observed mass; a planet's true mass exceeds its ob- 
served mass by at least one order of magnitude only if 
the planet's orbit is within sa 5.8° of being face-on. The 
currently observed orbital eccentricities of all the plan- 
ets studied here vary from nearly circular (w 0.01) to 
si gnificantly eccentri c (~ .53) . 

iLibert fe Henrarcfl (|2007l ) analyzed the proximity of 
particular multi-planet systems to strong mean motion 
resonance and found that HD 38529, HD 168443, HD 
190360 are dominated by non-resonant motions. For HD 
11964, the semimajor axis ratio of both planets typically 
exceeds 10, so a mean motion resonance is unlikely to 
dominate the evolution of that system. Of the 5 sys- 
tems studied, only HD 108874 is near a mean motion 
resonance of order < 10. 

3.2. Non-Keplerian Effects 

Because some of the above planets harbor semima- 
jor axes of just a few tenths of an AU, we must con- 
sider the possible effects on the dyn amics due to the 
pericenter precession cau sed by GR (|Ford et al.l 120001 : 
Adams & Laughlin 2006a, b). This precession is the most 
significant deviation from a Keplerian orbit at these or- 
bital periods and we use th e standard approximation 
(|Fabrvckv fc Tremaindl2007l) : 
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where n refers to the mean orbital motion, c the speed of 
light, M* the mass of the star, and where the subscripts 
"inner" and "outer" will henceforth refer to the inner 
and outer planets. In order to assess the possible effect of 
GR on our simulations, we can compare the value of wqr 
with the magnitude of the prec ession caused by se cular 
evolution with the outer planet (|Zhou fc Sunll200"3[ ): 
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and bs (ct) are Laplace coefficients, with a represent- 
ing the semimajor axis ratio such that a < 1. Equation 
^ is derived based on Laplace-Lagrange secular theory, 
accurate to second-order in the eccentricities. The plan- 
etary eccentricities in our simulations may attain values 
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close to unity, and higher-order expansions do not admit 
eigenvalue solutions such as those from Eq. ((2J (see, e.g., 
I Veras fe Armitagell2007D . Further, by definition, secular 
theory assumes constant semimajor axes. Hence, the es- 
timate of w sec for potentially near-resonant systems such 
as HD 108874 may be inaccurate. 

With the above caveats in mind, the ratio of these two 
precession rates, Xscc = w scc /wgr, provides a measure of 
the importance of GR in a system. Tabled] lists values of 
Xsec for representative values of orbital parameters and 
masses of all systems studied. Because Xscc ^> 1 only 
for HD 168443 and HD 108874, general relativity may 
play a significant role in the dynamics in HD 11964, HD 
38529 and HD 190360. Therefore, we explicitly add into 
Mercury the effect of GR according to the approxima- 
tion in Eq. ([TJ for calculating an additional set of 6500 
simulations for each of the three systems. 

For highly inclined systems, a planet may undergo cou- 
pled inclination and eccentricity oscillations through a 
secular exch ange of angular momen tum, known as Kozai 
oscillations. iKiseleva et al.l (|1998D gives the period for 
these oscillations - the "Kozai period" - as: 
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Because this period refers to the oscillations of the inner 
planet's argument of pericenter, the relative contribu- 
tion of the Kozai mechanism to GR can be expressed as: 
XKozai = WKozai/^GR = 2ir / P Koz uj gr . Table Q] reports 
representative values of XKozai for each system studied 
here. The values demonstrate that the Kozai mechanism 
could be significant (if the relative inclination exceeds 
w 40°) for HD 168443 and HD 108874, the two systems 
where secular perturbations from the outer planet dom- 
inates the effect of GR. 

Other possibly significant effects for Hot Jupiters are 
tidal deformation and stellar oblateness. In order to as- 
sess the magnitude of these additiona l effects, we de- 
rive ra tios from the expressions given bv lJordan fe Bakos 
([20M) : 
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where P* refers to the radius of the star and J2 is 
a oblateness constant present in the typical expres- 
sion for gravitational p otential (e.g., Eq. 6.218 of 
iMurrav fe Dermottlll999h . Both ratios Xtido and Xobi are 
only weakly dependent on the eccentricity of the inner 
planet. The ratios can be estimated only in a rough 
sense because of the generally unknown stellar Ji for 
the former and because of the unknown radii of several 
planets and the strong dependence on Ph m cr for the lat- 
ter. Further, for Xtidc, not shown are the dependencies 



on the internal structures of both the star and planet, 
which are unknown. Table Q] lists estimates for the ini- 
tial values of these ratio s by adopting J 2 = 2 x 10~ 7 
(|Pireaux fc Ro zclot 2003|) and assuming Pinner ~ Pjup- 
In almost all cases, Xtidc << 1 and Xobi << 1, so GR 
precession will dominate over tides and the effects from 
stellar oblateness. The only possible exception is for HD 
190360, with Xtide = 0.21. Hence, we neglect these latter 
two effects in our simulations. 

3.3. System Stability 

For all systems except HD 168443, over 95% of the 
nearly coplanar systems, in both the prograde and ret- 
rograde sense, remained stable over 1 Myr. Largely in- 
dependent of each planet's line-of-sight inclination, all 
systems demonstrate that as the relative inclination be- 
tween both planets approaches 90°, the systems become 
increasingly likely to undergo instability. FigureQ]graph- 
ically depicts the fraction of systems which remain stable 
over 1 Myr for all five systems. This figure broadly quan- 
tifies the effect that GR has on stabilizing systems with 
a close (a < 0.1 AU) planet. 

When Kozai oscillations are present, the maximum 
eccentricity achieved by the inner planet (ei nn cr,max) is 
lower in our numerical simulations which include GR 
than in the simulations which do not include GR. The 
Kozai mechanism induces angular momentum transfer 
which causes the inner planet's eccentricity, inclination 
and argument of pericenter to evolve. Because these 
effects are coupled, when GR causes a change in the 
rate of argument of pericenter, the effect affects the 
planet's eccentricity a s well. Both iLin et al.l (|2000h and 
iFabryckv fc Tremaind (j2007) relate the Kozai timescale, 
general rclativistic precession, and ei nner max through: 
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assuming that the inner planet is on an initially circular 
orbit and where the "0" subscript will henceforth refer 
to an initial value. We have generalized this equation to 
include an initial non-zero eccentricity (see the Appendix 
and Eqs. I6ll7l) . 

Equation ([5]) demonstrates that for an initial relative 
inclination high enough to trigger the Kozai mechanism, 
the maximum eccentricity of the inner planet decreases 
as u>gr increases. Hence, a high value of uigr provides a 
stabilizing mechanism for a system. 

We test the quality of the analytical approximation af- 
forded by Eq. (|6]) by using our simulation output. Apply- 
ing the values in Table [T]to Eq. ([7]) illustrates that for all 
systems except HD 190360, the inner planet's maximum 
eccentricity should not be reduced by GR and Kozai ef- 
fects. Fig. [5] plots this eccentricity difference vs. the 
initial relative inclination for simulation data (dots) and 
Eq. © (solid line) for HD 11964. The curve in Fig. [2] is 
based on a single representative value of XKozai; in real- 
ity, each dot corresponds to a different curve. The dots 
are generally below by the curve given in the analytical 
model, as expected. 

The relative contribution of Kozai-like oscillations can 
be identified by the libration of the argument of pericen- 
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ter of the inner planet (wi nner )- Asymmetric libration of 
this angle, particularly about the values 90° and 270°, 
is characteristic of the Kozai mechanism. Therefore, we 
compute the fraction of inner planet argument of peri- 
centers measured with respect to (as opposed to projected 
on) the invariable plane, for all systems. 

HD 108874 and HD 168443 perhaps best demonstrate 
the interplay between these angles because for those sys- 
tems (one near-resonant and one far from strong MMR) , 
XKozai ^ 1- Figures HE] illustrate the fraction of systems 
where Aw' librates about 0° (blue/dotted lines with 
triangles) or 180° (red/dashed line with squares), and 
where lui„ librates around 90° or 270° (black/solid line 
with dots). Prograde and retrograde coplanar systems 
predominately demonstrate antialignment and alignment 
of Aw' , respectively, for HD 108874, and alignment only 
in HD 168843. The solid lines peak as i rc i approaches 
90°. Therefore, the Kozai mechanism has a significant ef- 
fect at large relative inclination values, as demonstrated 
by the asymmetric libration of ui,. This libration would 
be masked by considering libration of Aw' alone, which 
would instead suggest that about half of the systems 
show alignment and the other half show antialignment 
at high relative inclinations. 

3.4. HD 11964 

Butl er et al.l (|2006l ) announced the existence of the 
multiple planet system around HD 11964. Soon after, 
iGregorvl (|2007bl) performed a Bayesian analysis of the 
same data and presented evidence for a third planet (of 
~ 0.21Afj) in the system at ~ 1.1 AU. This planet's 
existence would imply an unusually small jitter of less 
than 0.9 m/s. Further, it would require that the com- 
mensurability with a one year period is not the result 
of aliases with the annual obs erving pattern and/or im- 
precise barycentric correction ((Balucv 2008]). Therefore, 
there is only strong evidence for two planets given the 
current data, and we consider the dynamics of two-planet 
solutions. 

We ran 6500 N-body integrations of HD 11964 includ- 
ing GR and 6500 integrations not including GR. We 
present the simulation output in Tables [3] and [TUJ Table |3] 
indicates: 1) The median amplitude of the outer planet's 
current eccentricity does not exceed 0.075 for any incli- 
nation. The value of k c (= e c _ rn i n /e c . max ), indicates that 
the eccentricity variation achieved throughout the simu- 
lation is modest; k c > 0.83 in all cases with GR. Regard- 
less of the system's initial orientation, the planet's orbit 
does not become significantly more circular during its 
secular evolution. 2) In the coplanar prograde and retro- 
grade (general relativistic) cases, the eccentricity evolu- 
tion is quite modest (e.g. n b w 0.92(0.95), 0.89(0.93) and 
k c k, 0.92(0.94), 0.93(0.93)). 3) We find no indication of 
instability except when 60° < i re i,o < 120°. When i re i 
approaches 90°, the eccentricity of the inner planet be- 
gins to vary drastically. Incorporating GR precession de- 
creases the inner planet's eccentricity variation enough to 
nearly double the fraction of stable systems when 60° < 
*rci,o < 75° and 105° < i re i,o < 120°, and quintuple the 
fraction of stable systems when 75° < i re l,o < 105°. 

Table QJJ indicates: 4) In the coplanar prograde and 
retrograde cases, Aw' > 90°, indicative of large am- 
plitude libration (regardless of the inclusion of GR). 5) 
About three-quarters of systems with 45° < i re i,o < 75° 



and 105° < i ro i,o < 135° demonstrate asymmetric libra- 
tion of ujb, indicating that the Kozai mechanism does 
play a significant role at these high relative inclinations. 
6) Secular perturbations lead to complex dynamics at 
both low and high inclinations. For moderate relative 
inclinations (30° < i re i, < 45°, 120° < i re i,o < 135°), 
Aw' « 64° and is just as likely to librate about 0°, 180°, 
or some other value. These relative inclinations produce 
large RMS amplitudes, and in the prograde case, demon- 
strate preferential alignment of Aw' . 

3.5. HD 38529 

The planet HD 38529 b was announced in lFischer et al.1 
(l2001f) and the o u ter pl anet (HD 38529 c) was confirmed 
in iFischer et al.l (|2003D . HD 38529 c has a minimum 
mass within tenths of a Jupiter mass of the oft-cited 13 
-Wjup cutoff. In addition to orbiting planets, the host 
star exhibits IR excess, i ndicating the presence of a cool 
disk at wide separations (jMoro-Martm et "aT1l2007t l. 

TableslSlandllTlsummarize the results of all 13,000 sim- 
ulations for HD 38529. The first two tables demonstrate: 
1) For relative inclinations within 60° of coplanarity, the 
vast majority of systems are stable. 2) Almost all systems 
that are within 15° of i Te \,Q = 90° become unstable when 
neglecting GR; including GR increases the number of sta- 
ble systems by over a factor of 50 to 11% — 14%. 3) The 
outer planet's eccentricity remains relatively constant (to 
within 10%) in all GR-based cases (and all non-GR-base 
cases except for i re io values closest to 90°). Recall that 
GR primarily contributes to the pericenter change of the 
inner planet. Hence, the effect on the outer planet is in- 
direct, yet non-negligible at high relative inclinations. 4) 
The eccentricity variation of the inner planet increases 
gradually as i re i,o approaches 90°. Including GR slightly 
decreases the variability of the inner planet's eccentricity. 
5) For highly inclined, stable systems, the relative incli- 
nation can vary from less than 90° (prograde) to greater 
than 90° (retrograde). 

Regarding libration properties, Table [TT1 illustrates : 6) 
The RMS deviations of Aw' are high enough to be con- 
sistent with circulation in all relative inclination bins, 7) 
Assuming libration of Aw' , alignment is strongly pre- 
ferred in the coplanar prograde and retrograde cases. 8) 
GR's primary effect on the libration of Aw' is to cause 
preferential alignment of the angle at moderate relative 
inclinations (45° < i re i,o < 75° and 105° < i re i l0 < 135°). 
9) For these values of z ro i,o, asymmetric libration of Aw' 
is dominant, indicating the role of the Kozai mechanism. 

3.6. HD 108874 

HD 108874's two known planets (|Vogt et al.l 120021 . 
|2005[) have a orbital period ratio close to four, sug- 
gesting that they may be in or near a 4:1 MMR 
(|Gozdziewski et al.l [2001 iLibert fc Henrardl 120071 ). Our 
joint statistical and dynamical analysis is well suited to 
testing this possibility. First, we performed a standard 
MCMC analysis assuming that each planet traveled on 
an independent Keplerian orbit. As for the previous 
planetary systems, we used the standard MCMC output 
to generate ensembles of initial conditions for N-body 
integrations to test for dynamical stability and secular 
evolution of the eccentricity and inclinations. In princi- 
ple, these results might be affected by the planet-planet 
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interactions, particularly if the two planets are in a mean- 
motion resonance. Therefore, we performed a second 
analysis using full N-body integrations for both the sta- 
tistical and dynamical analyses. The prior probability 
distributions were similar to those in our previous anal- 
ysis with non-interacting planets. In other words, the 
priors are based on the period and velocity amplitude 
inferred from the Jacobi orbital elements. We adopt a 
stellar mass of 0.95M Q . Since the relative orientation of 
the orbits matters for the N-body simulations, we must 
add priors for the inclination of the orbital plane and 
longitude of nodes (where the orbit punctures the sky 
plane). We consider both: 1) a prograde coplanar con- 
figuration with an isotropic prior for the inclination of 
the system (i.e., a common inclination and a common 
longitude of ascending node for each system), or 2) an 
isotropic prior for the inclination of each orbit and a uni- 
form prior for the line of nodes. As before, we select 
subsets of the output based on the relative inclination 
to be used for long-term N-body simulations to inves- 
tigate stability and secular eccentricity and inclination 
evolution. The likelihood is evaluated based on compar- 
ing observations and the simulated stellar velocities (in 
a barycentric frame) evaluated at the time of each radial 
velocity observation. 

We again employ Markov chain Monte Carlo (MCMC) 
to sample from the posterior probability distribution for 
the masses and orbital parameters. The convergence rate 
of MCMC is often sensitive to the algorithm for propos- 
ing the next state (e.g., candidate transition probability 
distribution function). When performing Bayesian anal- 
ysis with N-body interactions, we replace the Gaussian 
random walk Metropolis algorithm adapted for analyzing 
radial velocity data sets (Ford 2006) with a differential 
evolution Markov chain Monte Carlo algorithm (DEM- 
CMC; ter Braak 2006). In this algorithm, each state of 
the Markov chain consists of an ensemble of initial con- 
ditions (known as a generation) . For our analysis of HD 
108874, each generation consisted of an ensemble of 1,024 
sets of initial conditions. For the initial generation, the 
initial conditions were drawn from a posterior sample cal- 
culated assuming independent Keplerian orbits (i.e., us- 
ing the standard MCMC approach; Ford 2006). For each 
of 10 5 subsequent generations, 1,024 new sets of initial 
conditions are proposed, based on combinations of multi- 
ple sets of initial conditions from the previous generation. 
The acceptance probability is chosen so as to converge to 
the desired posterior probability distribution. As before, 
we perform long-term N-body integrations for a subset 
of the initial conditions drawn from the latter portions 
of the Markov chain. 

The results of the DEMCMC+N-body analysis (in- 
cluding planet interactions) and the standard MCMC 
analysis (assuming independent Keplerian orbits) are ex- 
tremely similar (Table [5] For example, for coplanar, pro- 
grade systems, the MCMC (DEMCMC+N-body) analy- 
sis estimates the fraction of stable systems to be 97.6% 
(99.0%). Similarly, the statistics describing the extend of 
eccentricity evolution are Kb — 0.29 (0.29) and k c = 0.25 
(0.26). Thus, we conclude that our results for HD 108874 
are not sensitive to the effects of short-term planet-planet 
interactions or the assumption of independent Keplerian 
orbits. The exception is when sini becomes small. How- 
ever, if sini is small enough to cause significant short- 



term interactions, then the system is very likely to be 
rejected based the test for orbital stability. 

In order to explore the possibility of the system being 
in resonance, we perform a systematic search for the pos- 
sible libration of all resonant angles associated with the 
4:1, 7:2, 11:3, and 15:4 MMRs. These commensurabili- 
ties were chosen based on their proximity to an orbital 
period ratio of four. Each resonant angle has the form 

4> = .7lAout+i2Ai n +j3tI7 ou t+j4ro in +j5n ou t-|-j6^in, (9) 

where A, represents the mean longitude, w represents the 
longitude of pericenter, and ji are constants. We sam- 
ple four MMRs, up to order 11, in case any high- (> 10) 
order resonan ces may play a role in the s ystem's dynami- 
cal evolution ()Michtchenko et aLll2006bl ) . For each of the 
four MMRs, all combinations of ji are considered such 
that J2ji = and (js + je) is even. We performed this 
search for all 6,500 systems simulated for HD 108874. We 
flagged any angles which librated for at least to = 10 5 
yr. We discovered that less than 0.5% of all simulations 
studied exhibited libration, and even then only for a few 
10 5 yr. There is no apparent correlation with resonance 
and stability nor resonance and initial relative inclina- 
tion. Thus, HD 108874 is very unlikely to be in resonance 
based on the current RV data. 

Regardless of the possible influence of resonant or near- 
resonant terms, the characteristics of the secular evolu- 
tion (Table [6]) are striking: 1) Despite a median initial 
planet b eccentricity of 0.12, Kb < 0.30 for every rel- 
ative inclination bin, and approaches zero as irc^o ap- 
proaches 90°, 2) For ifg^o = 40° — 75°, the eccentricity 
of the planetary orbits approaches zero, 3) Unlike for 
the inner planet's orbit, the outer planet's orbit demon- 
strates an increasingly restrictive eccentricity range in 
the retrograde case as the planets' orbits become copla- 
nar. The values in Table [T^] reveal other interesting as- 
pects of the motion: 4) Aw' always demonstrates sym- 
metric libration, regardless of i re i,o- 5) Prograde coplanar 
systems are almost entirely anti-aligned, whereas retro- 
grade coplanar systems are entirely aligned. This shift 
from aligned to anti-aligned could indicate a significant 
contribution from near-resonant terms even though the 
planets are not strictly locked in a MMR. 6) Prograde 
coplanar systems feature the smallest libration ampli- 
tudes. 7) As demonstrated in Fig. [3j the system evolu- 
tion exhibits libration of ujf, in a robust Kozai-like fashion 
for 40° < « rc i,o < 140° except for the unstable systems 
prevalent in the 90° < i rc i.o < 120° regime. 

3.7. HD 168443 

Aft er the discov e ry of HD 168443 b (jMarcv et al.1 
I2001 . lUdrv et al.1 (12001 found a comp an ion b rown 
dwarf in the system! Weras fc Armitagg 's (|2007h nu- 
merical integrations of the system, assuming coplanarity, 
demonstrate that both planets undergo significant eccen- 
tricity evolution. They found that the two most apparent 
modulation frequencies are ~ 10 3 yr and ~ 10 4 yr, and 
that the apsidal angle clearly circulates. 

Our 6,500 simulations of this system are summarized 
in Tables 171 and IT3l Table [7] demonstrates: 1) Almost no 
stable systems can exist for 60° < £ re l,o < 120° and few 
systems remain stable when this range is extended by 15° 
in each direction. 2) At least 20% of all systems become 
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unstable for every relative inclination bin except for the 
perfectly coplanar case, where 13.0% of the systems sam- 
pled become unstable. 3) For the vast majority of stable 
systems, planet c's eccentricity variation is independent 
of the relative inclination. 4) For highly inclined, stable 
systems, the relative inclination can vary from less than 
90° (prograde) to greater than 90° (retrograde). Table HUl 
shows that: 5) For all relative inclinations, the libration 
amplitude of Aw' exceeds 90°, 6) Aw' is nearly always 
preferentially aligned and 7) Wf, exhibits strong Kozai- 
like behavior for i IC \ q values just under 60° and just over 
120°. 

3.8. HD 190360 

lUdrv et all (1200311 a nnou nced a planet in HP 190360 
and IVogt et al.l ( 20051 ) and iGozdziewski fe Migaszewskl 
( 2006) followed up with a detection and confirmation of 
an additional planet whose mass is approximately one 
Neptune-mass. Primarily because of the inner planet's 
small semimajor axis, GR has the greatest effect on this 
system (see Table [TJ. 

Tables El [9] and Q3] provide the summary output from 
the system's 6,500 GR-based and 6,500 non-GR-based 
simulations. Table [8] demonstrates: 1) An appreciable 
fraction of stable systems exist for all relative inclina- 
tions. However, for 60° < iro^o < 120°, ps 89% of systems 
without GR are unstable. Including GR, however, sta- 
bilizes nearly every one of those systems. 2) The outer 
planet's eccentricity varies by less than 1% in the vast 
majority of individual stable systems. 3) The presence 
of GR reduces the variation of the inner planet's eccen- 
tricity by up to a factor of 4. 

Table Q3] suggests: 4) Without GR, the lowest libra- 
tion amplitudes and the highest amplitude standard de- 
viations of Aw 1 occur when the relative inclination is 
15° — 30° offset from either the prograde or retrograde 
planarity. Including GR results in different conclusions; 
the system does not demonstrate small libration ampli- 
tudes nor asymmetric libration of Aw' . 5) In the copla- 
nar and retrograde planar cases, alignment of Aw' is 
preferred over anti-alignment by a factor of nearly four. 
6) With the presence of GR, the libration amplitudes 
all exceed 90°, commensurate with circulation. 7) The 
Kozai signatures are weaker than in any other system, 
corroborating the trend exhibited by the values of XKozai 
from Table [U 

3.9. HD 12661 

IVeras fc Fordl (|2009D performed a similar MCMC- 
based statistical analysis for HD 12661, but with broader 
relative inclination bins. For completeness, we com- 
pare some of their results with the outcomes presented 
here. HD 12661 is a two planet system where the inner 
(outer) planet properties are: M ps 2.3Afj up (1.8Mj up ), 
a w 0.8AU(2.4AU), and e ps 0.35(0.05). The authors 
demonstrated that the outer planet spends over ~ 96% 
of the time following an orbit that is more eccentric than 
the one currently observed. In the coplanar prograde 
case, all systems were stable, and in the coplanar ret- 
rograde case, almost all systems were stable. 16.4% of 
systems remained stable for 60° < i ro i,o < 90°, whereas 
only 0.2% did so for 90° < i ro i, < 120°. This asymme- 
try about 90° is significant and much greater than the 
slight asymmetries about 90° seen for HD 11964, HD 



38529, and HD 190360, but comparable to those in HD 
108874. Further, the alignment of Aw' for HD 12661 
shows a strong dependence on i re i,oj only ps 33% of copla- 
nar prograde systems show alignment of Aw', whereas 
fa 99.4% of coplanar retrograde systems show alignment. 
The only system studied here with a similar trend is HD 
108874, whose Aw' is nearly entirely antialigned in the 
prograde coplanar case, and entirely aligned in the ret- 
rograde coplanar case. 

4. DISCUSSION 

Hierarchical systems can be highly eccentric or in- 
clined (or even retrograde). Such large eccentricities 
and inclinations imply that numerical integrations are 
often necessary to model the dynamics of extrasolar 
systems. Analyti cal treatments of two-planet evolu- 
tion abound (e.g. iMurrav fc Dermottl [19991 iMorbidellil 
l2002t and references therein), and often rely on ap- 
proximations to the gravitational potential, known as 
"disturbing function s," The Laplace expa nsion of the 
disturbing function (jEllis fc Murray! I2000D relies on a 
Taylor expansion about zero eccentricities and inclina- 
tions and is limited in eccentricity by the Sundman 
criter ion (|Ferraz-M cllo 1993; ISidlichovskv fc Nesvornvl 
1994). Although effective for small Solar System bod- 
ies, this disturbing function mus t be used w i th care 
with respec t to extrasolar planets dV eras 2007]). iBeaugel 
(|1996ft and iBeauge fc Michtchenkol (|2003l ) developed a 
high eccentricity version of a disturbing function in the 
planar case which can reliably model the evolution of 
massive extrasolar planets with eccentricities up to pa 
0.5. Because quadrupole and octopole expansions are 
in terms of the semimajor axis ratio, their accuracy is 
diminished for weakly hierarchical systems. For non- 
copl anar systems, one can use sem i -numerical averag- 
ing (IMichtchenko fc Malhotral 12004 IMichtchenko et al.l 
2006br iMigaszewski fc G ozdzicwski 2009) or adaptions 
of Gauss' method (jTouma et al.ll2009t ). which are valid 
for high e and /. However, these studies neglect reso- 
nant and short-term perturbations. We find that these 
can be significant for some planetary systems, even if not 
in MMR (e.g. HD 108874). 

We can compare our secular e yolution results with 
those f rom other investigators. iBarnes fc Greenbergl 
(20061)) performed a dynamical study of 2-planet systems 
which did not incorporate uncertainties in the initial con- 
ditions and assumed that each system was coplanar and 
edge-on. With those assumptions, they found that sev- 
eral of these systems had at least one eccentricity that 
periodically obtains a small value. They claimed that 
an unexpectedly high frequency of systems with "near 
separatrix motion" implied that this property couldn't 
be due to planet scattering from initially circular orbits. 
Instead, we find that once we 1) focus our attention on 
two-planet systems with published RV data and signifi- 
cant secular evolution, 2) avoid systems near MMRs, and 
3) account for the uncertainties in the current orbital el- 
ements, only a modest fraction show significant secular 
eccentricity evolution with one planet returning to a near 
circular orbit for all viable inclinations. This fraction is 
equal to 2/6500 for HD 11964, 1/6500 for HD 38529, 
757/6500 for HD 108874, 27/6500 for HD 168443, and 
0/6500 for HD 190360, assuming that the planet which 
returns to or stays on a near circular orbit achieves an 
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eccentricity 1) range of at least 0.1, and 2) between 0.01 
and 0.05 at least once every 5 x 10 4 yrs. For HD 108874, 
if we specify that the eccentricity range must be at least 
0.2, then the fraction becomes 488/6500. 

In particular, Fig. [5] plots the value of k for the five 
systems analyzed here as a function of initial i re l,o- The 
figure demonstrates that the eccentricity variation of the 
inner planet typically increases as the initial i re l,o ap- 
proaches 9 0°, and that k can take on values from ~ 
0.10 - 0.90. iBarnes fc Greenber g (2006b) use a measure, 

e = 2min(v / u 2 + v 2 )/ (u max - u min + f max - Umin), where 
u = einnere utor sin (Aw) and v = e inner e ou ter cos (Am) to 
compare the minimum eccentricity value of either planet 
attained during evolution to the range of variation, as- 
suming the planetary orbits are coplanar. Their value of 
e for the near-MMR system HD 108874 (0.198) is com- 
parable to the value of k for near-coplanar prograde and 
retrograde values (see green/dashed line in Figure [5J. 
However, their e values for HD 38529 (0.44), HD 168443 
(0.219) and HD 190360 (0.38) are much less than the val- 
ues of k for the near-planar regimes. Thus, we find that 
hierarchical systems typically have more moderate eccen- 
tricity evolution than suggested bv lBarnes &: Greenbergj 
(2006b). Hence, we do not find a significant discrepancy 
between the secular evolution of hierarchical two-planet 
systems and simulations of planet scattering. We find 
that the eccentricities of the inner planets in these sys- 
tems vary by a factor of 2 or more for initial i vc \ o values 
cl oser to 90° than to 0° or 180°. 

iLibert fc Henrardl (|2006h have studied the secular evo- 
lution of HD 38529 and HD 168443. The secular eccen- 
tricity evolution of the outer massive planet in HD 38529 
is poorl y modeled by low-order L aplace-Lagrange secular 
theory ()Veras fc Armita gc 2007); N-body simulations in 
the coplanar, edge-on case indicate that the inner eccen- 
tricity evolves on a times cale of ~ 10 5 yr, and tha t the 
apsidal angle circulates. ILibert fc Henrardl ()2006f ) esti- 
mate that n b = 0.87(0.90) and k c = 0.999(0.999) for 
i los = 90° by using their high-order expansion (second- 
order Laplace-Lagrange theory in parenthesis). These 
values compare favorably with the values in the copla- 
nar row of Table [H with the difference attributed to 
our isotropic distribution of line-of-sight inclinations and 
choice of initial values for the orbital elements. We 
find slightly larger ecce ntricity evolution b ut qu alita- 
tively similar behavior. ILibert fc Henrardl (|2006f ) esti- 
mate that for ilos = 90° and using their high-order 
expansion (Laplace-Lagrange theory), Kb — 0.85(0.89), 
and k c — 0.83(0.91). The value of k c differs significantly 
from our estimate (0.49) in Table [7] We believe this 
difference is due to the the increased masses resulting 
from our isotropic line-of-slight inclination distribution, 
resulting in a larger extent of eccentricity evolution of 
the outer planet. 

Several researchers have considered the potential for 
additional small planets to exist in these systems, 
partic ularly those plane ts potentially in the habitable 
zone. iJones et al.l (|2006| ) considered likely habitable zone 
ranges, h z and habitability prospects (yes/no/maybe) for 
all 6 of the systems studied here. They found: HD 11964 
(h z = 1.6 - 3.1 AU, yes); HD 38529 (h z = 2.4 - 4.8 AU, 
no); HD 108874 (h z = 1.0 - 2.0 AU, no); HD 168443 
(h z = 1.1 - 2.3 AU, no); HD 12661 (h z = 1.0 - 2.0 AU, 



no); H P 190360 (h z = 1.0 - 2.0 AU, maybe). lErdi et al.l 
(2004) agree that planets in the habitable zone of HD 
38529 would likely have chaotic orbits. Similarly, both 
lErdi et al.l (|2004D and IBarnes fc Ravmondl ([20041 ) agree 
that the prospects for sta ble planets in the habi table zone 
of HD 168443 are slim. ISchwarz etail (|2007l) raise the 
possibility that dynamically stable Trojan planets of HD 
108874 b may exist in the habitable zone. 

In order to help address the question of whether HD 
190360 may admit habitable terrestrial planets, we per- 
form 2 additional sets of 1 Myr simulations with a disk of 
test particles (a good representation of terrestrial planets 
in systems with giant planets). We start from the out- 
put of MCMC simulations and randomly chose 2 sets of 
stable coplanar initial orbital parameters for the system. 
For each system, we superimposed a coplanar disk of 40 
terrestrial planets distributed in semimajor axis accord- 
ing to a Keplerian (power law of —3/2) distribution. We 
assigned the eccentricities of the planetesimals random 
values under 0.001 and randomized their mean anomalies 
and longitude of pcriastrons. The extent of the planetes- 
imal disk was [1.027, 2.022] A U (the same repres entative 
habitable zone values used bv I Jones" eFaDHool) for one 
system, and [0.128 x (1 + 0.01), 3.92 x (1 - 0.36)] AU 
(the representative range from the apocenter of the in- 
ner planet to the pericenter of the outer planet) for the 
other system. In the first (second) system, we find that 
91.4% (92.7%) of all terrestrial planets survive within 
the initial planetesimal disk. The surviving planets have 
modest eccentricities: the median and standard devia- 
tion of these values are 0.18 ± 0.12 (0.15 ± 0.14). Most 
initial instances of instability occur within 10 4 -10 5 yr. 
These results indicate that that HD 190360 can likely 
host a stable terrestrial planet in the habitable zone. 

5. CONCLUSION 

Due to limitations of real astronomical observations, 
often a significant range of planetary masses is consistent 
with observations. Therefore, it is necessary to investi- 
gate the dynamics of ensembles of planetary orbits and 
masses to accurately model the orbital evolution of exo- 
plancts. Hierarchical multi-planet systems demonstrate 
a wide variety of dynamical behaviors depending on the 
"i-LOS and *rci values, which are only weakly constrained 
by observations. Inclusion of GR in simulations of multi- 
planet systems with a Hot Jupiter may crucially affect 
the long-term stability, extent of eccentricity variation, 
and apsidal configuration directly. The eccentricity and 
inclination evolution of stable highly inclined systems are 
often dominated by Kozai-like oscillations, but can be 
limited by precession due to other planets or GR. 
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APPENDIX 

Because most of our simulations treat planets with non-zero eccentricities, we derive a more general relation than 
that of Eq. (0). The two conserved quantities for a two-planet system including GR based on a quadrupole-order 
potential are 2 : 

F = ~ Zl „2 =7 + ^2 {-2 - 3e 2 mcr + (3 + 12e 2 mcr - 15e 2 nncr cos 2 w ilmor ) sin 2 J re i >0 } , (1) 



e 2 

inner/ 



H 2 =z 3 (l- e 2 nncr ) cos 2 J re i,o, (2) 



where 



3G 2 Af*M inncr (M* + M innor ) 

Zl = 12 3 . ( 3 ) 

"inncr L 

(j-M*-M; nner -Mouter a innei / a\ 

Z2 ~ Ml + M~ 8fl 3 (1 _ e 2 n3/2 ' <■ ^ 
""■outer V 1 ^outcr/ 1 



Noting that Z\/z2 — 6Pkoz^scc = 127r/xKozai7 one can then compute F' = Fjz-x and H' 2 = H 2 / ' z% from the initial 
conditions and then solve the following implicit equation for ei nncr . max : 

2 ■ 3/ 2 \ ^ ^ ^ e inncr,max / ' ' \ ( 27T . 

cos z reL0 = - (1 - e inncrjmax ) - -- 2 — 2 - I — ) (6) 

° ° e inner,max e inner,0 \ . / 1 — e 2 n \1~ e- I V AKo 



inner, y inner, max , 

Because Eq. ([5]) may admit no solutions for certain values of / r ei,0i one must consider the physically plausible solutions. 
When the inner planet's initial eccentricity equals zero, no Kozai librations occur for XKozai < 27r/3. More generally, 
as a function of ei nncrj o, no Kozai librations occur when: 

XKozai < * V , ' ■ (7) 

^e 2 . /1 -e 2 



inner, y inner, 



Equation J7]) demonstrates that the critical XKozai for which Kozai oscillations will occur is a weak function of ei nncri o 
so that such oscillations should always occur for XKozai > 2-7T when < ei, m er.o J$ 0.79. The equations, however, are 
just approximations, because they assume 1) secular evolution only (no change in semimajor axis), 2) a Hamiltonian 
truncated to the quadruple-order term, 3) a doubly-averaged Hamiltonian which eliminates short-period terms, and 
4) a fixed outer orbit whose orbital plane does not vary. 



Note that Eq. (17) of lFabrvckv fc Tremaine] <2007l) should read U H' 2 = (l - e? 



MCMC Catalog 



TABLE 1 
Ratio of Effects 



System Name 


Xsec 


XKozai 


Xtide 


Xobl 


HD 11964 


1.02 


6.43 


0.027 


0.0036 


HD 38529 


0.79 


6.11 


0.017 


0.006 


HD 168443 


94.4 


628.7 


0.0003 


0.0013 


HD 190360 


0.14 


1.05 


0.21 


0.0022 


HD 108874 


1869.7 


11068.1 


0.000003 


0.00029 



Note. — Values of the ratios, Xscc = ^scc/^ghi 
Xtidc = ^tido/^GR and Xobl = ^obl/^GR- The table 
values indicate that general relativity may have an ap- 
preciable effect on the dynamical evolution of HD 11964, 
HD 38529 and HD 190360, and that in most cases for 
each of the 5 systems, effects from stellar oblateness and 
tidal effects are weak compared to those from general rel- 
ativity. These values are computed from representative 
values of the minimum masses and orbital parameters of 
each system. 



TABLE 2 

HD 11964 WITHOUT GENERAL RELATIVITY 





Stable 




e;, curve 


Kb 


^■c,med 


e c curve 




^r,med 




0° 


100.0% 


0.28 


0.028,0.15,0.30,0.45,0.57 


0.92 


0.068 


[0.012,0.055,0.11,0.17,0.21] 


0.92 




0° 


- 15° 


100.0% 


0.30 


0.032,0.16,0.32,0.48,0.61 


0.90 


0.069 


[0.012,0.055,0.11,0.17,0.21] 


0.92 


10.5° 


15° 


- 30° 


100.0% 


0.31 


0.035,0.18,0.35,0.53,0.67 


0.75 


0.071 


[0.012,0.052,0.11,0.16,0.21] 


0.92 


23.9° 


30° 


- 45° 


100.0% 


0.30 


0.035,0.18,0.36,0.55,0.69 


0.47 


0.074 


[0.012,0.052,0.11,0.16,0.20] 


0.94 


38.7° 


45° 


- 60° 


100.0% 


0.30 


0.045,0.22,0.44,0.66,0.83 


0.32 


0.068 


[0.012,0.052,0.11,0.16,0.21] 


0.91 


53.1° 


60° 


- 75° 


48.8% 


0.28 


0.045,0.23,0.47,0.70,0.89 


0.24 


0.074 


[0.012,0.052,0.11,0.16,0.20] 


0.86 


64.4° 


75° 


-90° 


4.4% 


0.23 


0.048,0.24,0.48,0.72,0.92 


0.17 


0.038 


[0.0083, 0.045, 0.092, 0.14, 0.18] 


0.57 


83.7° 


90° 


- 105° 


4.6% 


0.30 


0.078,0.27,0.50,0.74,0.93 


0.19 


0.053 


[0.0083, 0.038, 0.075, 0.12, 0.21] 


0.67 


95.8° 


105° 


- 120° 


40.2% 


0.29 


0.045,0.23,0.46,0.69,0.87 


0.25 


0.067 


[0.0083, 0.048, 0.098, 0.15, 0.24] 


0.87 


116.6° 


120° 


- 135° 


99.6% 


0.30 


0.045,0.23,0.46,0.69,0.87 


0.31 


0.067 


[0.012, 0.055, 0.108, 0.17, 0.21] 


0.90 


127.1° 


135° 


- 150° 


100.0% 


0.30 


0.038,0.20,0.40,0.59,0.75 


0.45 


0.071 


[0.012,0.052,0.11,0.16,0.21] 


0.93 


141.9° 


150° 


- 165° 


100.0% 


0.29 


0.032,0.17,0.33,0.49,0.63 


0.75 


0.066 


[0.0083, 0.045, 0.092, 0.14, 0.18] 


0.91 


156.1° 


165° 


- 180° 


100.0% 


0.28 


0.035,0.17,0.35,0.52,0.66 


0.89 


0.076 


[0.012,0.055,0.11,0.16,0.21] 


0.93 


169.4° 



Note. — Summary of stability, eccentricity evolution, and relative inclination evolution for HD 11964 without the 
inclusion of GR arranged by bins of relative inclination (Column 1). Column 2 displays the percent of stable systems 
out of 500, Columns 3, 6 and 9 display median starting values, and Columns 4, and 7 display the 5, 25, 50, 75 and 95 
percentiles of the time-averaged eccentricity across all stable integrations within the given initial relative inclinations. 
Columns 5 and 8 display the mean ratios of the minimum vs. maximum eccentricities obtained throughout the 
simulations. 



TABLE 3 

HD 11964 WITH GENERAL RELATIVITY 





'rel 


Stable 


&b,med 


ej, curve 








e c curve 






^r,rned 




0° 


100.0% 


0.28 


0.028,0.15,0.30,0.44,0.57 


0.95 


0.068 




0.012,0.052,0.11,0.16,0.21 




0.94 




0° 


- 15° 


100.0% 


0.28 


0.028,0.15,0.30,0.45,0.57 


0.93 


0.068 




0.012,0.052,0.11,0.16,0.21 




0.93 


10.6° 


15° 


- 30° 


100.0% 


0.28 


0.028,0.15,0.30,0.45,0.57 


0.83 


0.068 




0.018,0.098,0.20,0.29,0.37 




0.94 


23.3° 


30° 


-45° 


100.0% 


0.28 


0.035,0.18,0.36,0.54,0.69 


0.61 


0.068 




0.012,0.055,0.11,0.16,0.21 




0.94 


38.5° 


45° 


-60° 


100.0% 


0.28 


0.042,0.21,0.41,0.62,0.79 


0.36 


0.068 




0.012,0.055,0.11,0.16,0.21 




0.90 


52.4° 


60° 


- 75° 


88.2% 


0.28 


0.045,0.23,0.46,0.68,0.87 


0.27 


0.066 




0.012,0.052,0.11,0.16,0.21 




0.86 


66.9° 


75° 


- 90° 


27.8% 


0.34 


0.048,0.24,0.48,0.72,0.92 


0.28 


0.064 


[0.0083, 0.048, 0.098, 0.15, 0.21] 


0.83 


81.8° 


90° 


- 105° 


23.6% 


0.33 


0.052,0.24,0.48,0.72,0.92 


0.28 


0.068 


0.0083, 0.048, 0.098, 0.15, 0.19 


0.84 


96.7° 


105° 


- 120° 


83.6% 


0.27 


0.045,0.23,0.45,0.67,0.85 


0.26 


0.068 


[0.012,0.055,0.11,0.16,0.21 




0.84 


113.3° 


120° 


- 135° 


100.0% 


0.28 


0.042,0.21,0.41,0.62,0.78 


0.34 


0.068 


[0.012,0.055,0.11,0.17,0.21] 


0.90 


127.3° 


135° 


- 150° 


100.0% 


0.28 


0.042,0.20,0.41,0.61,0.77 


0.59 


0.068 


[0.012,0.052,0.11,0.16,0.20] 


0.94 


141.5° 


150° 


- 165° 


100.0% 


0.28 


0.038,0.19,0.38,0.57,0.73 


0.83 


0.068 


[0.012,0.058,0.12,0.18,0.22] 


0.94 


156.8° 


165° 


- 180° 


100.0% 


0.28 


0.028,0.15,0.30,0.45,0.57 


0.93 


0.068 


[0.012,0.052,0.11,0.16,0.20] 


0.93 


169.1° 



Note. — 



Same as Table [2] but incorporating the effects of GR. 
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TABLE 4 
HD 38529 WITHOUT GENERAL 
RELATIVITY 



Table 4 is located on www.dimitrivcras.com 



TABLE 5 
HD 38529 WITH GENERAL 
RELATIVITY 



Table 5 is located on www.dimitrivcras.com 



TABLE 6 
HD 108874 



Table 6 is located on www.dimitrivcras.com 



TABLE 7 
HD 168443 



Tabic 7 is located on www.dimitrivcras.com 



TABLE 8 
HD 190360 WITHOUT GENERAL 
RELATIVITY 



Table 8 is located on www.dimitrivcras.com 



TABLE 9 
HD 190360 WITH GENERAL 
RELATIVITY 



Table 9 is located on www.dimitrivcras.com 
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TABLE 10 

HD 11964 (WITH GENERAL RELATIVITY) 



*rel 



aligned 



antialigned 



other 



RMS Aro' 



0° 
15° 
30° 
45° 
60° 
75° 
90° 
105° 
120° 
135° 
150° 
165° 



0° 

- 15° 
-30° 

- 45° 
-60° 

- 75° 

- 90° 

- 105° 

- 120° 

- 135° 

- 150° 

- 165° 

- 180° 



68.6% 
65.8% 
60.6% 
39.6% 
37.6% 
50.0% 
50.0% 
60.9% 
42.3% 
31.3% 
34.2% 
60.8% 
57.0% 



(68.4%) 
(66.8%) 
(66.2%) 
(57.0%) 
(51.2%) 
(54.2%) 
(54.7%) 
(40.7%) 
(40.0%) 
(35.0%) 
(55.0%) 
(61.4%) 
(78.8%) 



59.8% 
60.6%' 
54.4%' 
41.2%' 
19.0%' 
15.6% : 
36.4% : 
13.0% 
17.3%' 
19.2%' 
34.0%' 
48.6%' 
52.0%' 



26.8% (30 
31.0% (32 
34.4% (32 
21.6% (39 
30.4% (30 
34.0% (35 
27.3% (41 
21.7% (55 
37.8% (47 
34.9% (44 
25.4% (39 
32.8% (36 
34.6% (20 



8%) [6.2%] 
4%) [8.8%] 
6%) [10.2%] 
0%) [13.2%] 
4%) [6.6%] 
.4%) [9.0%] 
7%) [13.6%] 
9%) [21.7%] 
4%) [8.4%] 
.6%) [8.0%] 
8%) [15.2%] 
6%) [12.6%] 
6%) [11.2%] 



4.6% (0.8%) 34.0% 
3.2% (0.8%) 30.6% 
5.0% (1.2%) 35.4% : 
38.8% (4.0%) [45.6%] 
32.0% (18.4%) [74.4%] 
16.0% (10.4%) [75.4%] 
22.7% (3.6%) [50.0%] 
17.4% (3.4%) [65.2%] 
19.9% (12.7%) [74.3%] 
33.7% (20.4%) [72.7%] 
40.4% (5.2%) [50.8%] 
6.4% (2.0%) [38.8%] 
8.4% (0.6%) [36.8%] 



91.4° ± 5.0° 
91.8° ±4.3° 
91.3° ±6.5° 
63.6° ± 25.4° 
64.3° ± 25.0° 
74.4° ±21.5° 
81.2° ± 12.6° 
75.4° ±20.2° 
75.1° ±21.5° 
62.7° ±27.4° 
64.6° ± 25.4° 
92.4° ± 5.8° 
93.0° ±5.6° 



(95.7° ±4.9°) 
(94.9° ±4.9°) 
(93.8° ±4.4°) 
(91.5° ±8.9°) 
" ±20.2°) 
± 16.6°) 
± 11.4°) 
±11.8°) 
±17.1°) 
±21.7°) 
(91.1° ±9.8°) 
(94.4° ±4.4°) 
(95.0° ±4.8°) 



(78.8 
(80.3 
(84.0 
(83.7 
(78.5 
(77.7 



95.2° 
96.3° 
96.6° 
70.9° 
65.4° 
66.7° 
75.7° 
70.0° 
69.8° 
64.9° 
[69.9° 
[95.8° 
[96.5° 



±5.5° 
±5.4°' 
±7.1°' 
± 24.2° 
±26.7 
±26.4 
± 19.9 
±22.2 
±25.8 
±26.9 
± 24.6° 
±7.8°] 
±7.1°] 



Note. — Summary of libration of Aw' about aligned configurations (Column 2), antialigned configurations (Column 3) and 
asymmetric configurations (Column 4) arranged by bins of relative inclination (Column 1) for HD 11964. Column 5 displays 
the root mean squared (RMS) libration and its standard deviation across different initial conditions about the most prevalent 
libration center. Values in parenthesis are computed when GR is included in the simulations, and values in square brackets 
refer to the libration of cji n measured with respect to the invariable plane. 



TABLE 11 

HD 38529 WITH GENERAL RELATIVITY 



Table 11 is located on www.dimitrivcras.com 



TABLE 12 
HD 108874 



Table 12 is located on www.dimitrivcras.com 



TABLE 13 
HD 168443 



Tabic 13 is located on www.dimitrivcras.com 



TABLE 14 
HD 190360 



Table 14 is located on www.dimitrivcras.com 
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Fig. 1— The percent of stable systems for HD 11964, HD 38529, HD 108874, HD 168443 and HD 190360 as a function of relative 
inclination when the effect of general relativity is included (blue/dashed lines and squares) and when it is not (black/solid lines and dots). 
The data points correspond to 15° -wide relative inclination bins along with a coplanar prograde bin. 
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Initial Relative Inclination (Deg) 

Fig. 2. — The maximum eccentricity achieved by the inner planet of HD 11964 for stable systems as a function of initial relative inclination 
(dots). Overlayed is the solid curve predicted from theory (Eq. |6j which predicts the maximum eccentricity attained by planets affected 
by general relativity, assuming an initial representative inner planet eccentricity of 0.28. Dots appear above the curve at low inclinations 
because in this regime, the eccentricity oscillations are not dominated by Kozai effects, but rather by classical Laplace-Lagrange secular 
theory. 
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Fig. 3. — The percent of simulated HD 108874 stable systems whose apsidal angle projected on the invariable plane (Aro') is predominantly 
aligned (blue/dotted line with triangles) and antialigncd (red/dashed line with squares). Overlayed is the percent of systems whose inner 
planet argument of pericenter measured with respect to the invariable plane librates around 90° or 270° in a Kozai-likc fashion (black/solid 
line with dots). The data points correspond to 15°-wide relative inclination bins along with a coplanar prograde bin. Bins with no symbols 
or lines contain no stable systems. 
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Fig. 4. — Same as Fig. [3]but for HD 168843, a system with two planets not near MMR. 




Fig. 5. — The averaged eccentricity variation (ft = e m i n /e m ax) of the inner planet as a function of initial relative inclination for HD 
11964 (black/solid lines with triangles), HD 38529 (blue/dotted lines with squares), HD 108874 (green/dashed lines with dots), HD 168443 
(rcd/dot-dashed lines with diamonds) and HD 190360 (magenta/triple dot-dashed lines with asterisks). The data points correspond to 
15°-wide relative inclination bins along with a coplanar prograde bin. Bins with no symbols or lines contain no stable systems. 



